AD-A248  845 


E^ni  >ci.WKMr  kkA»&I^ICATlON  |J 

Jnclasslfled 


ICURITY  classification  a 


ECLASSIFICATION !  DOWNC 


(FORMUKC  OAGANIZATiON  re^t  nuvbe 


ITATION  PAGE 


form  Approvttf 
OMf  No.  07M-Of  M 


•  ME  OF  fCRFORMING  ORGANIZATION 
resttr  Polytechnic  Institute 


>ORES$  (Oty,  SUU,  and  ZIP  Coda) 
100  Institute  Road 
Worcester,  HA  01609 


AVE  OF  FUNDING /SPONSORING 
RSANIZATION 

!:  Group 


DORESS  roty.  Sutt,  and  ZIP  Coda) 

>,  Army  Materials  Technology  Laboratory 
ertown,  HA  02172 


TlF  (includa  Sacurity  Oasuftcation) 

Jltrascund  NDE  of.  Adhesive  Bond  Integrity:  A  Quantitative  Measure 


ERSONAL  AUTHORS) . 

-udwlg,  Relnhold  and  Sullivan,  John  M. 


8&  OFFICE  Symbol 
(It  apfilicabla) 

SLCMT-MRM 


tb  RESTRICTIVE  MARKINGS 

None  _ _ _ 


1.  DISTRIBUTION /AVAILABILITY  OF  REPORT 

Unlialted 


S  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


7«  NAME  OF  MONITORING  ORGANIZATION 

13-S.  Army  Materials  Technology  Laboratory 


7b  ADDRESS  (Oty,  ttata,  and  ZIP  Coda) 
Watertovm,  MA  02172 
« 


9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

DAAL04-90-C-0024 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 

PROJEa 

TASK 

ELEMENT  NO. 

NO. 

NO. 

n/PE  OF  REPORT 
Final  Report 


13b  TIME  COVERED 
FROM  8/17/90  TO 8/1 8/91 


14  DATE  OF  REPORT  (Yaar,  Idonttt.  Day)  TS  PAGE  COUNT 

1992,02.28  I  81 


COSAT  I  CODES 


GROUP  SUB-GROUP 


IB  SUBiECT  TERMS  (Continue  on  ravana  it  nacasury  and  idantify  by  block  numbar) 
Ultrasonic  NDE,  Contact  Transducer,  Signal  Processing, 
Chirp-Z  Transform,  Adaptive  Filtering,  Least  Squares, 
Deconvolution,  Bondline  Transfer  Function 


ISTRaCT  (Continue  on  reverie  it  nacatsary  and  idantify  by  block  numbar) 


ee  attached  sheet. 


1)2  4  16  Ot>i 


92-09853 


iitiiiiiiiiimn 


ilSTRlBUTION /availability  OF  ABSTRAH  121  ABSTRAH  SECURITY  CLASSIFICATION 

UNCLASSIFIED/UNLIMITED  □  SAME  AS  RPT  □  DTiC  USERS  I  Unclassified 


NAM'  OF  responsible  INDIVIDUAL  I22b  TELEPHONE  (inc/uON  Aree  Code/  22c  OFFICE  SYMBOl 

Robert  Anastasl  _  _ 1  (617)  923-5241  SLO!T-MR>( 


brm  1473,  JUN  S4  Pnaioutadrtionaaraobaolata 


nSjffwwPTBHfynfswiflfflBCfff 


I:  SUMMARY 


The  objective  of  this  research  was  to  investigate  and  develop  a  coupled 
approach  (analytical,  numerical  and  experimental)  to  the  ultrasonic 
nondestructive  evaluation  of  adhesive  bond  integrity.  Results  of  tiiese 
studies  were  directed  toward  nondestructive  evaluation  (NDE)  of  die  integ¬ 
rity  of  adhesive  bonds  and  bondlines  for  advanced  composites  and  multi¬ 
layered  materials  by  the  U.S.  Army.  Such  studies  are  needed  in  both  the 
manufacturing  phases  of  products  in  order  to  assure  quality  as  well  as 
during  the  operating  lifetime  of  the  products  in  order  to  predict  water 
infiltration  or  bond  deterioration  and  hence  prevent  failure. 

The  developed  multiple-staged  model  started  from  a  mathematical 
description  of  the  ultrasound  propagation  through  an  inhomogeneous, 
isotropic  or  anisotropic  solid  with  appropriate  boundary  conditions  for  the 
transmitter/receiver  unit  coupled  to  the  material  under  test.  Analytical 
solutions  were  employed  to  initially  test  and  calibrate  the  numerical 
fcimulations.  This  numerical  approach  was  configured  to  be  flexible  and 
realistic  enough  to  investigate  a  wide  variety  of  bond  configurations  on  the 
computer. 

Our  strategy  was  targeted  at  die  coupling  of  these  numerical  simula¬ 
tions  of  the  underlying  physical  processes  with  the  experimental  data 
gathered  in-house  at  MTL.  The  interaction  of  transducer  signals  with 
different  simulated  bondline  configurations  such  as  disbonds,  pure  bonds, 
and  bonds  of  varying  thickness  were  examined.  The  resulting  synthetic 
data  was  compared  to  experimental  measurements.  The  correlations  of  the 
experimental  and  numerical  trials  were  extended  to  include  feature  extrac¬ 
tion  capability  related  to  interfacial  bond  thickness.  Therein,  individual 

characteristic  features  lepreseoting  the  bond  thickness  were  numerically  _ 

isolated  and  overlayed  with  the  experimentally  observed  signals.  ?.r 

This  novel  multi-stage  approach  addressed  the  generation  and  propaga-  * 

tion  of  elastic  wave;>,  the  wave  interaction  at  the  interface  and  advanced  |  ^ 

signal  processing  of  experimental  data.  |— _ 1  ***^^** 

An  executive  summary  addressing  each  component  of  the  original  !  ■> _  _ _ 

proposal  is  followed  by  detai  led  chapters  of  each  work  item.  :  “  1 ‘  r  i  ^  i  •» 
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I:  EXECUTlVliSUMxMARY 


PROJECT  DESCRIPTION 
A.)  Background 

The  integrity  of  bonded  structures  is  of  paramount  importance  in  the 
safe  and  reliable  operation  of  milifaiy  equipment.  The  Patriot  rocket 
requires  adhesive  attachment  of  the  ceramic-dome  to  the  Kevlar  ring. 
Similarly,  titanium  components  must  be  bonded  to  phenolic-based 
materials.  Helicopter  rotor  blades  are  multilayered  composites  bonded 
together  Cobra  deck  panels,  missle  radomes,  mines,  projectiles,  and 
explosive  cartridges  arc  samples  ot  advanced  military  equipment  requiring 
adhesive  bonds  and  bondlines.  The  operational  readiness  and  security  of 
these  units  depend  to  a  large  extent  on  the  integrity  of  the  interfacial  bonds. 

Nondestructive  evaluation  of  the  bonds  and  bondlines  has  experienced 
verj  limited  success  [1].  NDE  methods  are  needed  to  determine  the 
strength  of  the  bonds  in  situ.  Unfortunately,  no  NDE  method  has  demon- 
stnited  the  ability  to  quantitatively  state  the  strength  of  a  bond.  Several 
NDE  techniques  have  been  applied  to  access  adhesive  bond  quality.  They 
include  ultrasonics,  acoustic  emissions,  radiography,  holography,  nuclear 
magnetic  resonance,  eddy  current,  and  thermal  imaging. [2]  Of  these  ac¬ 
cepted  NDE  methods  only  ultrasonics  appears  to  retain  a  reasonable  proba¬ 
bility  of  success  in  the  boiidline  application  [3].  Thermal  imaging  requires 
the  .x)mposite  or  radome  shell  to  be  of  exbemely  small  thickness  Other¬ 
wise,  tlie  thermal  image  of  tlie  interface  is  completely  masked  by  the  diffu¬ 
sion  of  heat  tiirough  the  substrate.  Nuclear  magnetic  resonance  has  shown 
Some  Mccess  in  the  detection  rf  water  vithin  bonds  but  it  has  not  been 
ablt  0  detect  disbc»nds,  foreign  matte»-  ur  weak  bonds.  Ultrasonic  pulse- 
echo  techniques  have  been  shown  to  reliably  detect  total  disbond  regions. 
Some  investigators  indicate  that  the  strength  of  the  adhesive  layer  can  be 
correlated  to  the  attenuation  coefficient  and  velocity  of  sound  in  the  mate¬ 
rial  [4].  Interfacial  and  horizontally  polarized  shear  wave  measurements 
have  demonstrated  the  ability  to  discriminate  bond  strength  [5].  However, 
the  test  configurations  were  highly  restrictive  and  the  practical  application 
of  the  technique  is  undefined  as  yet.  Leaky  Lamb  waves  technique  uses  an 
oblique  incident  wave.  Usually  the  waves  reflect  off  the  interface,  howev- 
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er,  certain  frequencies  excite  plate  waves  in  the  structure.  These  leaky  waves 
have  been  shown  to  interfere  with  the  expected  reflection  patterns  [6].  Deter¬ 
mining  the  correct  frequency  for  destructive  interference  patterns  is  not  reli¬ 
able.  Researchers  have  swept  the  frequencies  to  find  a  zone  of  interference. 

If  one  examines  multiple  frequency  responses  of  the  same  test  configuration 
wnall  subtleties  of  the  bond  integrity  can  be  inferred.  This  ability  to  detect 
weak  bonds  is  still  marginal.  However,  it  does  indicate  that  if  one  could  cor¬ 
relate  the  multiple  values  obtained  from  the  complex  input  signals  then  the 
probability  of  characterizing  the  system  configuration  would  increase. 

B.)  Research  Accomplishments 

The  basis  of  the  research  plan  was  the  development  of  a  multiple  attack 
strategy  as  shown  in  Figure  1.  Building  on  our  past  experience  into  the  analyt¬ 
ical  study  of  transient  elastic  waves  radiating  into  an  elastic  half-space  [8],  a 
precise  model  was  developed  which  served  as  a  calibration  for  our  numerical 
formulation.  This  numerical  model  solves  the  transient  elastodynamic  equa¬ 
tion  of  motion  subject  to  realistic  boundary  and  initial  conditions  in  two-di¬ 
mensional  (x,y)  space  and  in  three-dimensional  (axisymetric)  space. 

Depending  on  the  employed  probe,  a  longitudinal  or  shear  wave  contact 
transducer  have  been  simulated.  The  boundary  and  initial  conditions  can  be 
set  over  a  wide  range  of  practically  relevant  bounds.  Furthermore,  since  the 
model  discretizes  the  general  stress  equation  of  motion,  anisotropic  and  inho¬ 
mogeneous  material  parameters  can  be  taken  into  account.  It  is  this  feature  of 
matenal  inhomogeneiiy  that  makes  our  numerical  modeling  approach  particu¬ 
larly  suitable  for  bondiine  inspection.  A  realistic  transducer  response  signal 
was  incorporated  into  the  numerical  model  and  tested  against  the  analytical 
theoiy,  die  propagation  of  tlic  acoustic  pulse  was  monitored  throughout  the 
lest  specimen.  )Mth  the  displacement  field  given  at  discrete  instances  in  time 
it  was  possible  to  “freeze”  the  acoustic  pulse  prior  to  reaching  the  bonded 
area.  The  bond  line  itself  can  be  studied  in  a  broad  fashion  depending  on  geo¬ 
metric,  material,  and  density  parameters.  Figure  2  shows  some  of  the  bond 
conngurations  tested  using  our  numerical  system  and  experimentally  tested  at 
MTL.  The  simulated  signals  are  currently  being  stored  and  compared  with 
practical  measurements  to  train  a  neural  net  The  initial  feature  targeted  for 
extraction  from  the  data  is  bond  thickness  which  is  the  main  focus  in  the  sec¬ 
ond  phase  of  this  project  (MTL  #  DODAJ.S.  Army  -  DAAL04-91-C-(X)54). 
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Figure  1  •  Multipte  Attack  Strategy  for  the  Soludoo  of  Adhesive  Bond  Integrity 


Figure  2a  -  Double  area  discontinuity  model  -  Uniformly 
distributed  microscopic  disbonds  resulting  from  surface 
contamination  -  Taken  from  Thompson. 


Figure  2b  -  Captured  displacement  field  in  upper  substrate 
is  applied  to  numerous  bondltne  configurations 

The  use  of  our  numerical  model  has  allowed  us: 

a)  to  gain  crucial  insight  into  the  physical  processes  of  ultrasound/ 
bondline  interaction.  This  view  of  the  interior  region  and  the  physics 
of  the  system  are  evident  in  Figure  3.  Note  how  the  numerical  model 
shows  the  longitudinal,  shear,  Rayleigh,  and  head  waves  propagating 
and  reflecting  within  the  material.  This  transparant  window  is  un¬ 
available  in  an  experimental  investigation. 
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b)  to  generate  synthetic  data  which  can  be  correlated  with  practical 
measurements  in  order  to  interpiete  field  data  of  relevance  to  die  U.S. 
Army; 

c)  although  not  implemented  completely  in  this  phase  I  study,  the  nu¬ 
merical  model  will  allow  us  to  optimize  the  NDE  inspection  process 
by  determing  optimal  transducer  configurations,  aperture  sizes,  fre¬ 
quency  settings,  etc. 


U.  lurf**  plot  for  plane  «iiin  with.  1/4  inch  r-0.8  iis 


Our  numerical  apprcacli  was  closely  coupled  with  tlie  in-house 
experimental  investigations  at  MTI^.  The  simulated  transient  outputs 
(A -scans)  were  directly  compared  to  the  MTI.  measurements.  The 
result  section  of  this  report  shows  the  excellent  agreement  between 
the  two  investigation  modes.  Additionally,  the  numerical  approach 
aided  in  the  development  of  new  or  modified  signal  processing  meth¬ 
ods  such  as  the  Chirp-Z  approach  currently  under  investigation  at 
MTL  [12]. 
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C.)  Woiic  Statement  of  Research 

WPI  in  conjunction  with  MTL  explored  the  ultrasound  bondline  NDE  integ¬ 
rity  from  all  three  accepted  modes  of  scientific  endeavor  -  i.e.  experimental, 
analytical,  and  numerical  modes. 

•  Experimental  Investigations 

a. )  MTL  testing 

The  dominant  focus  was  in  experimental  measurement  of  samples  at 
MTL.  These  bond  samples  underwent  extensive  analysis  using  longitudi¬ 
nal  and  horizontally  polarized  shear  waves.  The  Chirp-Z  transformation 
techniques  were  applied  to  enhance  the  obtained  measurements.  The 
experimental  output  was  saved  in  multiple  formats.  The  unfiltered  trans¬ 
ducer  signals  (A-scans)  as  well  as  hardware  Chirp-Z  transformed  outputs 
were  stored  for  all  test  configurations. 

b. )  Feature  Analysis 

The  rationale  for  this  procedure  was  that  a  specimen  prepared  with  a 
prescribed  flaw  or  geometry  exhibits  a  characteristic  signal  response  as 
given  by  such  fealiues  as  frequency  content,  maximum  amplitude,  rise 
time,  etc.  [13].  We  have  tested  experimentally  and  numerically  multiple 
layer  specimens  with  variety  of  thicknesses.  These  signatures  are  inputs 
to  a  neural  net  software  system  which  is  being  trained  to  output  discrete 
thickness  values  of  the  interface.  The  training  is  ongoing  with  encourag¬ 
ing  results  and  is  a  dominant  component  of  our  phase  n  effort.  The  neural 
net  program  was  obtained  by  MTL. 

•  Mathematical  Problem  Formulation 

Based  on  the  proposed  multiple  attack  strategy  the  analytical  formula¬ 
tions  of  the  problem  were  thoroughly  investigated  as  to  their  applicability 
for  realistically  modeling  th<",  physical  process.  Initial  and  boundary  con¬ 
ditions  were  examined  relative  to  their  influence  on  the  underlying  equa¬ 
tions.  Approximations  to  the  governing  equations  were  derived  and  inves¬ 
tigated. 

•  Numerical  Modeling 

A  finite  element  model  of  the  elastic  wave  equation  was  implemented 
on  a  mid-size  computer  (DecStation  5000  model  200).  The  code  includes 
pre-  and  post  processing  options  such  as  adaptive  mesh  generation  and 
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graphics  routines.  The  user  identifies  the  physical  properties  of  the  various 
layers  (density,  wave  velocities)  and  a  desired  resolution.  The  system  auto¬ 
matically  discretizes  a  domain  with  individualized  nodal  resolutions  in  each 
material.  The  successful  two-dimensional  modeling  (x,y)  was  extended  to 
account  for  full  three-dimensional  geometries  (axisymmetric). 

•  Verification 

The  verification  of  our  numerical  investigations  was  two  fold.  Using 
prescribed  test  configurations  our  simulations  were  compared  to  existing 
analytical  solutions  yielding  graphically  indistinguishable  results.  Equally 
important,  numerical  simulations  were  compared  to  experimental  field  tests 
gathered  at  MTL.  The  result  section  shows  that  the  numerical  simulations 
reproduced  fliat  observed  experimentally  with  fidelity. 


D.)  Summary 

This  woilc  presents  a  concise  overall  approach  to  the  analysis  of  bond-line 
mtegrity  as  measured  via  ultrasound  NDE.  Based  on  the  previously  outlined 
multiple  attack  strategy,  a  clear  mathematical  recipe  was  presented  that  de¬ 
scribes  the  primary  components  of  the  coupled  physical  phenomenon  of 
ultrasonic  waves  impinging  on  an  adhesive  interface.  It  is  expected  that  this 
novel  approach  will  increase  the  general  comprehension  of  ultrasonic  interfa¬ 
cial  interacUon  for  nondestructive  testing  purposes.  Moreover,  the  numerical 
analysis  approach  was  directly  tested  against  analytical  studies  given  certain 
simplifying  constraints  and  against  cxpoimental  results  gathered  at  MTL.  In 
each  situation  the  simulated  results  predicted  that  observed  in  the  analytical  or 
experimental  investigative  mode  with  fidelity. 
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2  ELASTODYNAMIC  WAVE  PROPAGATION 


2.1  Governing  Equations 


The  governing  equation  is  the  balance  of  linear  momentum  (Cauchy’s  first  law  of  motion) 
which  in  indicial  notation  can  be  expressed  as 


+  F. 


(2.1) 


where  is  the  stress  tensor,  F^  is  the  body  force  vector,  p  is  tiie  density  and  u^  is  the 

displacement  vector.  The  Einstein  summation  convention  is  assumed  where  /^1,2,3 
representing  the  three  coordinate  directions  x,  y  and  z  respectively.  (  ),j  represents  a 
spatial  derivative  with  respect  to  the  /  direction. 


In  order  to  solve  this  equation  for  displacements  it  is  necessary  to  express  the  stress 
tensor  in  terms  of  displacements.  A  Kelvin  model  constitutive  relation  is  chosen  to  allow 
viscous  damping.  [15] 


= 


(22) 


where  Cjjy  is  the  elastic  constitutive  matrix,  is  the  damping  constitutive  matrix, 

is  the  strain  matrix  and  is  the  rate  of  strain  tensor  A  single  dot  above  a  variable 

signifies  a  first  derivative  with  respect  to  time,  two  dots  a  second  derivative.  Choosing 
a  Maxwell  model,  or  a  more  complicated  model  for  that  matter,  makes  it  difficult  to  solve 
for  the  stress  directly.  It  is  assumed  that  the  strains  are  small  and  that  tiie  material 
behaves  in  a  linear  elastic  manner.  The  strain'displacement  and  rate  of  strain-rate  of 
displacement  relations  are  as  follows 


(23) 


The  nature  of  the  motion  can  be  categorized  into  three  types  of  waves;  longitudinal,  shear 
and  surface.  Longitudinal  (or  primary)  waves  are  characterized  by  particle  motion  that 
is  parallel  to  the  direction  of  propagation.  Shear  (also  called  transverse  or  secondary) 
waves  are  characterized  by  motion  that  is  perpendicular  to  the  direction  of  propagation. 
Surface  waves  have  motion  with  components  in  each  direction  but  whose  magnitude 
decreases  exponentially  away  from  the  surface.  Surface  waves  on  a  free  surface  are 
known  as  Rayleigh  waves.  Waves  at  solid-solid  boundaries  are  refened  to  as  Stonely 
waves. 
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Figure  2.1:  Wavefronts  generated  by  a  point  or  line  source. 


Figure  (2.1)  shows  the  wave  propagation  for  a  point  force  on  an  infinite  half-space.  All 
the  waves  propagate  perpendicular  to  Aeir  wavefronts.  The  schematic  only  identifies 
wavefront  positions:  it  does  not  reflect  the  wave  amplitudes  which  may  be  zero  at  some 
points  on  the  wavefront  The  longitudinal  wave  travels  with  greatest  speed  which  is 
expre-ssed  as'*  f-  fox  an  isotropic,  linear,  homogenous  medium  where,  iiandl/ 

are  Lamp's  constants  and  p  is  the  material  density.  The  shear  wave  speed  is  defined  as 

C,  =  The  Rayleigh  wave  has  a  velocity  between  86.2  ard  95.5  percent  of  the  C, 
depending  on  Poisson’s  ratio.  If  the  half-space  is  viewed  from  above,  the  Rayleigh  wave 
has  A  linear  wavefront  for  plane  strain  and  a  circular  wavefront  for  axisymmctric.  The 
head  wave  is  a  shear  wave  that  is  produced  by  the  longitudinal  wave  at  the  free  surface. 
The  mode  conversion  is  necessary  to  satisfy  the  stress  free  boundary  condition.[161 
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Transducer 


Figure  22:  Expected  wavefronts  for  finite  width  transducer. 


Figure  (2.2)  shows  the  expected  wavefronts  for  a  finite  width  transducer.  The  applied 
force  is  spread  uniformly  over  the  transducer  width.  The  discontinuity  at  the  edge  of  the 
applied  force  produces  wavefronts  similar  to  those  generated  by  a  point  force.  The  width 
of  the  transducer  only  contributes  to  the  flat  longitudinal  wave. 
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3  FINITE  ELEMENT  FORMULATION 


3.1  Numerical  Formulation 

Finite  element  solution  techniques  for  partial  differential  equations  arc  roost  commonly 
derived  by  variational  principles  or  by  the  method  of  weighted  residuals.  Each  method 
attempts  to  satisfy  the  differential  equation  by  approximating  the  solution  with  a  set  of 
known  functions  multiplied  by  a  set  of  arbitrary  constants. 

-  E 


Where  are  the  known  functions,  or  shape  functions,  and  the  unknown  coefRcients  are 
the  displacements  t/^.  The  capitalized  subscript  Af  refers  to  the  number  of  discrete 
points  p  used  to  approximate  the  solution  [18]. 

The  variational  approach  is  based  on  minimizing  a  scalar  functional  derived  from  the 
original  differential  equation.  When  solving  equilibrium  problems  of  elasticity  this 
functional  represents  the  potential  energy  of  the  system.  In  the  static  case,  the 
minimization  produces  an  extremum  of  the  functional  guaranteeing  convergence  of  the 
soliitio'i.  However,  applying  the  variational  method  to  the  elastodynamic  wave  equation 
prod,  .  cs  only  a  stationary  value  not  an  extremum  of  the  functional  [10],  therefore 
convc:gencc  is  not  guaranteed,  which  is  to  be  expected  considering  the  hyperbolic  nature 
of  the  wave  equation. 

The  weighted  residual  approach  will  be  presented  for  its  sinpler  mathematical  nature,  ie. 
knowledge  of  variational  calculus  is  not  required.  The  method  of  weighted  residuals  deals 
directly  with  the  differential  equations.  Substituting  the  assumed  solution  for  the  exact 
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solution  into  equation  (2.6),  the  differentia’  equation  will  usually  not  be  satisfied  over 
the  entire  domain.  The  resulting  error  is  known  as  the  residual. 


Rt  =  *  Pi  ~  pVat^M  *  0  (3^) 

The  spatial  derivative  does  not  appear  on  the  11^  ,  and  terms  because  diey  are 

^atial  constants.  This  error  is  then  forced  to  be  zero  on  average  over  the  domain  by 
setting  the  volume  integral  of  the  product  of  a  set  of  weighting  functions  Wj^  and  the 

residual  R^  equal  to  zero. 

*  D^U^yNy  *  F^-  pO^Ny]dV  =  0  (3J) 

V  -0  il 

This  formulation  produces  an  independent  equation  for  each  degree  of  freedom.  Although 
the  equations  can  be  solved  at  this  point,  Green’s  theorem  of  integration  by  parts  is  used 
on  the  first  two  terms  of  equation  (3.3).  This  practice  reduces  the  order  of  the  derivative 
on  the  assumed  solution,  allowing  simple  linear  shape  functions  to  be  used,  and  also 
expresses  the  boundary  conditions  explicitly. 

*  /,  *  /» 

*  fy^F.dy  -  fw^O^N„dy  =  0 


Note  that  the  exact  solution  has  been  reinserted  into  the  surface  integrals  and  that  is 
the  unit  outward  normal  on  the  surface.  By  combining  the  surface  integrals  and 
substituting  for  the  stress  tensor  using  equation  (2.2),  it  becomes  obvious  that  they 
represent  the  surface  tractions  at  the  boundary. 
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(3^) 

7|  is  the  surface  traction  vector  at  the  boundary.  Equation  (3.4)  can  be  further 
simplified  by  moving  die  constant  terms  out  of  die  integrals.  According  to  the 

Galerldn  mediod  of  weighted  residuals  setting  the  weighting  functions  equal  to  the  shape 
functions  yields 

*  lfN,D^N^dy]0^ 

*  t/ Wt 

Equation  (3.6)  is  applicable  to  the  entire  domain.  The  equation  is  effectively  discretized 
into  finite  elements  by  the  appropriate  choice  of  shape  functions  Nj^. 

3.2  Plain  Strain  and  Axisymmetric  Formulation 

In  order  to  avoid  the  complexities  of  a  full  three  dimensional  formulation,  the  sinqiler 
two  dimensional  plain  strain  and  axisymmetric  models  are  formulated.  Advantages 
include  reduced  total  degree  of  freedom  models,  simpler  mesh  generation,  and  ease  of 
graphical  presentation  of  results.  The  computational  strategics  developed  in  section  5  are 
not  limited  in  any  way  to  a  two  dimensional  formulation  and  can  be  easily  extended  to 
a  full  three  dimensional  formulation. 

The  plain  strain  formulation  models  an  object  infinitely  thick  in  the  z  coordinate  direction. 
This  eliminates  normal  strains  in  the  z  direction  as  well  as  shear  strains  in  the  x-z  and  y-z 
directions.  The  transducer  is  now  modeled  as  an  infinite  strip  and  any  voids  become 
infinite  tubes.  The  waves  generated  are  cylindrical  about  the  z  axis.  As  the  cylindrical 
waves  expand  the  amplitude  of  the  wave  decreases  because  its  energy  is  ({Head  over  a 
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larger  area.  This  phenomenon  is  referred  to  as  geometric  dispersion.  The  actual 
geometric  dispersion  for  a  circular  transducer  is  spherical  and  as  a  result  A-scans  from 
diis  formulation  do  not  compare  well  with  experimental  data  quantitatively. 

The  axisymmetric  formulation  models  a  body  of  revolution  under  axisymmetric  loading 
conditions.  This  formulation  models  the  physical  transducer  extremely  well  because  the 
geometric  dispersion  is  spherical.  The  formulation  does  not  work  well  for  voids  that  are 
off  of  the  centerline  which  are  modeled  as  circular  rings.  If  the  enphasis  of  the  analysis 
is  on  the  wave  interaction  between  layered  mediums  the  axisymmetric  formulation  is  the 
most  practical. 


The  fourth  order  tensors  of  equation  (3.6)  are  difficult  to  deal  with  and  can  be  reduced 
to  second  order  tensors  by  index  contraction  [19].  For  the  plain  strain  case  the  andD^ 
matrices  become 


An  axis  of  symmetry  may  be  chosen  in  an  or  y'  direction  if  the  stiffness  and  damping 
matrices  transform  such  that  =  0  and  =  0.  The  shape 

luijction  matrix  is  defined  as 


0  1^2  0  ...  0 

0  N,  0  AT,  ...  0 


(3.8) 


The  derivatives  of  the  shape  functions  become 

assuming  that  the  displacements  have  been  contracted  according  to 
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dx 

0 

A 

dN^ 

U 

if 

dN^ 

dy 

dx 

_  w,  .  awj 

-  0  -ST  0 

N  ay  ay 


[  dy  cbc  dr  dy  dx  j 

[t/,1  ... 

where  the  first  subscript  is  the  direction  and  the  second  is  the  node  number. 

Equation  (3.6)  can  now  be  cast  in  the  more  familiar  form 

\M](j  +  [D]U  +  [K]U  ^  F^*  F, 

with  the  following  definitions. 

(Ml  =  fN^pN^df  F,  =  fNJ’,didy 

A  A 

fDI  =  dxi,  F,  = 

[K]  -  fN^CiiSy  ibdy 


(3.10) 


(3.11) 


(3.12) 


where  [Af]  is  the  mass  matrix,  [D]  is  the  damping  matrix,  [JiT]  is  the  stiffness  matrix, 
are  the  nodal  body  forces  and  F,  arc  the  nodal  surface  forces.  The  volume  integrals 

have  becii  reduced  to  area  integrals  due  to  tlte  two  dimensional  formulation  in  which  a 
unit  depth  in  the  z  direction  has  been  assumed. 

For  the  axisymmetric  formulation,  is  the  same  and  C^,  and  take  the 
following  forms. 
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<^1122 

^1122 

0 

^1111 

^1122 

^1122 

0 

<^1122 

®2222 

^2233 

0 

*  Da  ® 

<^1122 

^ - - 

^«3QB  ZZS3 

0 

^1122 

^2233 

^2222 

0 

^1122 

^*233 

0 

0 

0 

0 

^12U 

0 

0 

0 

aiVj 

dNy 

0 

0  ... 

L 

0 

dz 

dN, 

dl 

av. 

dz 

dN, 

0 

0 

2 

0 

dr 

dr 

dr 

N.  » 

iV, 

A'i 

1 

0 

Z 

0  ... 

L 

0 

r 

r 

r 

^N^ 

dN, 

dN^ 

dN^ 

3r 

dz 

dr 

••• 

di 

a^ 

dz  , 

(3.13) 


(3.14) 


In  order  for  aii  axisymmetric  formulation  u  be  valid,  the  material  must  be  at  least 
transversely  isotrofic  with  the  material  axis  of  symmetry  parallel  to  the  z  axis  [1 IJ.  The 
expressions  presented  for  the  [C]  and  [f]  matrices  have  been  ^propriately  reduced. 

«  J InrNiF^drdz 

A 

F,  »  f2nrNJ,dS  (3.15) 

5 

The  integrations  in  equations  (3.i2)  and  (3.15)  can  be  calculated  numerically  using  Gauss 
quadrature.  The  shape  functions  for  common  element  types  and  the  procedure  for 
Gauss  quadrature  can  be  found  in  many  textbooks  [20,21]. 


Equation  (3.11)iS  now  defined  by 

[Af]  =  JlnrNj^pNydr^ 

A 

[D]  =  f  lurNj^D^N^fdrdz 
[A"]  -  j  2nrA(j  Cjfty  drdz 
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33  Time  Integration 


Equation  (3.11)  represents  a  system  of  linear  second  order  dMerendal  eqtiations  in  time. 
Numerical  solution  schemes  are  based  on  discretizing  time  into  intervals  and  satisfying 
the  system  only  at  those  time  slices.  Various  schemes  include  the  central  difference 
method,  the  Newmark  method,  the  Houbolt  mediod  and  the  Wilson-6  method.  The 
difference  between  these  schemes  is  determined  by  how  diey  estimate  Ae  displacements, 
velocities  and  accelerations  within  each  time  step  [22].  Although  these  methods  are 
commonly  referred  to  as  ffnite  difference  methods  in  time,  Zienldewicz  has  shown  that 
they  are  equivalent  to  a  finite  element  formulation  using  the  weighted  residual  method 
applied  to  an  interval  2Ar  if  the  shape  functions  are  chosen  appropriately  [23]. 

These  methods  fall  into  two  groups,  implicit  and  explicit  Explicit  implies  that  the 
solution  at  time  t  +  At  der^ends  only  on  known  values  at  previous  time  steps.  Implicit 
implies  that  the  solution  at  time  t  *  At  depends  on  known  values  at  previous  dme  steps 
as  well  as  values  at  the  current  time  step.  The  main  advantage  of  the  explicit  methods 
is  that  the  system  of  equations  is  decoupled.  The  disadvantage  is  that  small  time  steps 
arc  required  for  stability.  Implicit  schemes  are  stable  for  larger  time  steps  but  require  a 
matrix  to  be  either  inverted  or  triangularized  at  the  first  time  step.  Due  to  the  large  size 
of  the  matrices  involved,  the  explicit  central  difference  method  is  employed. 

Using  the  central  difference  method  the  acceleration  term  is  estimated  as 

U‘  =  -J— [C/'*Ar  .  217'  +  U'-^']  (3.16) 

(AO" 

assuming  that  At  remains  fixed.  The  velocity  term  was  estimated  using  a  backward 
difference  to  avoid  a  matrix  inversion  and  maintain  an  explicit  solution  [23] 

l/‘ =  J-[I7‘ -  £7'-^']  (3.17) 

At 
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Substituting  for  U‘  and  U*  using  equations  (3.16)  and  (3.17)  into  equation  (3.11)  and 
solving  for  U**^,  yields  Ae  following  fully  explicit  equation  neglecting  body  forces. 


«  Af[Ml-*(Af(F;  -  (iq£/')  +  {D](U*-^  -  I/')}  + 


(3.18) 


3.4  Mass  Matrix  Diagonalization 


The  mass  matrix  resulting  from  Ae  mtegration  defined  m  (3.12)  and  (3.15)  is  called  Ae 
consistent  mass  matrix.  It  can  be  evaluated  using  Gauss  quadrature  numerical  mtegration. 
However,  to  maintain  an  explicit  solution  Ae  mass  matrix  is  Aagonalized.  The  mass 
matrix  can  be  diagonalized  wiAout  a  significant  loss  of  accuracy  by  using  Ae  following 
formula  on  an  elemental  basis  [24].  . 


m; 


EE«; 


(3.19) 


where  Ae  superscripts  d  and  c  stand  for  Aagonal  and  consistent  mauices  respectively. 
Alternatively  Ae  Aagonalized  mass  matrix  can  be  evaluated  Arectly  by  moving  Ae  Gauss 
quatjratiue  points  to  Ac  node  points.  This  approach  works  well  wiA  Ae  plane  strain 
forimtJanon  but  fails  for  Ae  axisymmetric  case  when  Ae  elenjent  is  on  Ae  centerlme 
(r  -  0) .  After  diagonalization,  Ae  matrix  inversion  is  reduced  to  Averting  each  Aagonal 
component  mdividually.  If  a  central  Afference  formula  is  used  to  estimate  Ae  velocity. 

solvmg  for  would  require  Ae  sum  of  [M\  and  [D] ,  which  would  no  longer  be 

diagonal,  to  be  Averted.  The  [D]  matrix  can  not  be  Aagonalized  because  Ae  sum  of  Ae 
components  will  always  equal  zero. 
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3^  Stability  Constraints 


The  disadvantage  of  the  explicit  method  is  that  it  is  only  conditionally  stable.  The 
constraint  governing  the  stability  is  known  as  the  Courant  constraint 

Courant  =  q  =  C_ — i  1.0  (3.20) 

Where  is  the  longitudinal  wave  velocity  of  the  material  of  that  element  and 
Ax^  is  the  length  of  the  shoitest  edge  for  a  rectangular  four  noded  element  For  a  three 
noded  triangular  element  or  a  four  noded  parallelugram,  Ax^  will  be  the  shortest 

distance  from  any  node  to  a  non-adjacent  side.  Every  element  must  satisfy  equation 
(3.20)  for  the  solution  to  be  stable.  In  the  explicit  method  information  flows  across  one 
elemenf  per  time  step.  Therefore  the  minimum  rate  at  which  information  flows  is  the 
mininiuni  Ax  divided  by  Utc  time  step  At.  If  that  rate  is  slower  than  the  numerical 

solution  can  not  approximate  the  actual  wave  motion  and  the  solution  will  diverge,  which 
is  the  constraint  of  equation  (3.20).  [26] 

3.6  Initial  Conditions 

The  initial  conditions  needed  to  start  equation  (3.18)  are  the  full  displacement  vectors  at 
t  -  -At  and  r  =  0  The  first  displacements  calculated  are  for  t  =  Af .  Alternatively,  if  the 
acceleration,  velocity  and  displacemenl  vectors  are  known  at  r=0  a  Taylor  series 
expansion  can  be  used  to  approximate  the  displacement  at  f  =  At. 

=  t/o  +  Af(/®  +  |(Af)*t/®  (3.21) 

If  only  the  initial  velocity  and  displacement  vectors  are  known  a  backwards  Taylor  series 
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expansion  can  be  used  to  fold  the  displacements  at  t  =  -Af  as  follows  [11]. 

U~^  =  ~  (3.22) 

Equation  (3.18)  can  then  be  used  at  subsequent  time  steps.  Although  non-zero  initial 
conditions  are  acceptable,  it  is  assumed  that  the  medium  is  initially  at  rest  with  no 

internal  stresses.  The  initial  conditions  then  become  U~^  =  0  and  t/®  =0.  If  a  forcing 

function  is  applied  that  is  zero  at  r  »  0  the  solution  of  equation  (3. 18)  at  r  =  Ar  will  not 

produce  any  displacements  and  can  be  avoided  by  substituting  for  F'  in  equation 

(3. 18)  and  increasing  the  total  dmc  at  each  time  step  by  At  and  setting  the  response  atr  =  At 
equal  to  zero. 
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4  COMPUTATIONAL  CONSIDERATIONS 


Upon  first  inspection,  it  would  seem  that  implementing  the  finite  element  solution  could 
be  done  in  a  similar  manner  to  the  simple  static  case.  However,  it  becomes  apparent 
quickly  that  there  are  several  problems  The  extremely  short  wavelengths  of  the 
ultrasonic  waves  force  the  finite  element  mesh  to  become  extremely  dense.  Due  to  the 
large  number  of  degrees  of  freedom  of  the  dense  mesh,  the  memory  required  to  store  the 
matrices  is  excessive.  The  difficulty  is  compounded  by  the  requirement  of  a  time  step 
size  that  is  proportional  to  the  element  spacing  tinougb  the  Courant  constraint.  The 
number  of  time  steps  necessary  to  complete  an  analysis  quickly  grows  into  the  thousands. 
Corse<]ucntly  execution  times  Chi»  easily  stretch  into  months.  This  section  presents  a 
computer  implementation  strategy  that  svoids  these  pitfalls,  creating  a  useful  tool  for  the 
investigation  of  ultrasonic  phenomena. 

4.1  Memory  Storage 

Fo'  large  ineshcs  the  memory  required  to  store  the  stiffness  and  damping  matrices  quickly 
becomes  prohibitive  b&:ause  thei*  size  is  of  the  order  of  the  number  of  degrees  of 
l.;m  squared.  Foi  c.\a.-iiple,  <>'  the  mesh  has  1,000,000  degrees  of  freedom,  the  full 
sti^'fness  matrj<  v-d)l  ncr.c  4,000  gigabytes  of  memory  using  single  precision  For  a 
banded  matrix  solver,  assuming  a  500  by  1000  node  mesh,  the  memory  reduces  to  8 
gigabytes.  However  using  a  sparse  matrix  solver  only  requires  the  storage  of  the  non- 
zciO  values,  reducing  the  memory  to  72  megabytes  for  the  stiffness  matrix.  The  damping 
matrix  is  of  cqval  size  increasing  the  total  to  144  megabytes.  Although  tins  amount  could 
be  handled,  it  requires  the  epu  to  swap  memory  back  and  forth  from  secondary  storage 
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which  is  extremely  inefficient  and  time  consuming.  An  alternative  method  for  storing 
those  matrices  is  pursued. 


One  of  the  advantages  of  using  the  explicit  formulation  is  that  the  global  matrices  are  not 
assembled.  Bathe  and  Wilson  [22]  suggest  that  the  displacements  can  be  computed  using 
an  element  by  element  approach.  Displacements  are  multiplied  by  the  elemental  matrices 
and  the  contributions  are  assembled  in  the  new  displacement  vector.  Further  savings  can 
be  achieved  if  only  the  unique  elemental  matrices  are  saved.  Although  the  element  by 
elenvent  strategy  reduces  the  memory  required  to  the  absolute  minimum  it  creates 
redundant  operations.  Assuming  a  rectangular  linear  mesh,  each  node  will  be  visited  by 
four  elements.  There  are  eight  multiplications  for  each  degree  of  freedom  per  element 
visit  or  64  total  multiplications  per  node.  The  assembled  matrix  multiplication  is 
completed  by  eighteen  multiplications  per  degree  of  freedom  resulting  in  a  total  of  36 
multiplications  pei  node  For  uniformly  rr  .M.co  linear  triangular  elements,  there  arc  72 
multiplications  for  the  element  by  clement  approach  and  only  28  for  the  assembled 
matrix.  For  a  S-dimcnsional,  8  noded  brick  formulaticn  tiicrc  are  576  .lultiplications  for 
the  element  by  element  and  243  for  tlic  assembled  natrix. 


*  . 

■  •  UnlQM  El«n«n( 

.  -IMqivNodo 


Figure  4.1:  Plane  strain  unique  node  and 
eJement  locations. 


B  -  Uniqua  Bamant 
•  -UnlquaNoda 


Figure  4.2:  Axisymmetric  unique  node 
and  element  locations. 


This  redundancy  can  be  eliminated  by  identifying  and  storing  unique  rows  of  the 
assembled  matrices.  A  node  by  node  approach  can  then  be  utilized  to  reduce  the  number 


24 


uf  operations  to  the  minimum.  A  unique  node  location  is  determined  by  the  properties 
of  the  surrounding  elements  and  by  their  relative  positions.  For  simple  layered 
geometries,  with  homogeneous  material  properties  within  each  layer,  an  a  priori  method 
of  determining  the  unique  nodal  locations  is  shown  in  Figures  (4.1)  and  (4.2)  for  the  plain 
strain  and  axisymmetric  configurations,  respectively.  The  number  of  unique  nodes  for  the 
layered  model  can  be  calculated  using  the  following  formulas. 

Number  of  unique  nodes  =  6(NL)  +  3  (Plane  Strain) 

(4.1) 

Number  of  unique  nodes  =  2(NL  l)iNR)  (4xisymmetric) 

Where  NL  is  the  number  of  layers  and  NR  is  the  number  of  nodes  in  the  r  direction.  The 
plane  strain  formula  is  independent  of  die  number  of  nodes  in  the  x  or  y  directions.  The 
axisymmetric  formula  is  independent  of  the  number  of  nodes  in  the  z  direction  but 
dependent  on  the  number  of  nodes  in  the  r  direction  The  memory  required  to  store  the 
stiffness  and  damping  mairioes  for  a  5(Ki  by  1000  ncvle,  thice  layered  model  is  reduced 
Crom  144  Mbytes  to  12  kbytes  for  plane  strain  and  2,3  Mbytes  for  axisymmetric.  Each 
unique  node  is  given  an  integer  flag  to  associate  it  with  its  unique  row.  The  remaining 
non-unique  nodes  aie  assif.ned  the  .flag  number  for  the  appropriate  conesponding  unique 
node  Tiie  mass  matrix  can  be  stored  in  the  same  fashion  because  the  unique  nodes  in 
Figures  (4.1)  and  (4  2)  are  also  unique  mass  nodes. 

ibc  local  stiffness  matrices  arc  only  calculated  at  the  unique  element  locations  shown  in 
Figures  (4.1)  and  (4.2).  The  unique  clcnents  are  flagged  in  a  similar  manner  as  the 
unique  nodes.  Row  assemblage  is  accomplished  by  visiflng  each  unique  node,  calculating 
the  unique  element  flags  for  the  four  sunoundir.g  elements  and  assembling  the  appropriate 
values  from  the  local  stiffness  matrices  into  the  unique  rows,  one  for  each  degree  of 
freedom.  The  process  is  exactly  the  same  as  in  conventional  finite  element  methods 
except  that  onh  two  rows  of  the  assembled  matrix  are  saved. 

The  node  by  node  form  of  cquatiai  (3.18)  can  be  written  as  follows. 
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•c"-  -  «•)}  «  2l/i  -  C4‘'  («) 


The  bar  over  tfie  displacements  indicates  a  vector  quantity  and  (  indicates  the 
appropriate  row  of  the  matrix.  [  .  indicates  only  the  diagonal  conq>onent  of  the 

appropriate  row,  making  a  scalar  quantity. 

4.2  Solver  Efficiency 

The  bulk  of  the  time  used  by  the  solver  can  be  attributed  to  two  processes.  The  first  is 
the  actiud  solution  to  the  differential  equation  The  second  is  the  determination  of  the 
element  coimec  ivity  necessary  to  place  local  node  contributions  in  the  correct  global 
locations.  To  develop  the  most  efficient  solver  possible  both  of  these  processes  are 
streamlined. 

Equation  Reduction 

The  first  step  in  minimization  of  the  actual  solution  factors  equation  (4.2)  for  the 
displacements.  The  resulting  coefficients  are  constant  and  are  premultiplied  to  avoid 
repeating  the  same  multiplication>  every  time  step.  The  factored  equation  becomes 

iC*'  •  f,'  ♦  mL®*  ♦  ml  v"*' 

with  the  following  definitions. 

[JCC  “  *  2[/U 

ml  =  -  t/]„ 

where  [/]  is  the  identity  matrix.  If  damping  is  not  included,  equation  (4.3)  becomes 
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14“  -  F,'  *  mil/  - 14“ 


(«) 


with  the  following  definitions. 

f;  -  Ar»[Mii«,F; 

mL  =  -  2m^ 


(*J6) 


Note  that  the  U*~^  term  is  no  longer  a  vector  quantity  in  equation  (4.5). 

Automatic  Connectivity 

The  time  expended  on  the  element  connectivity  can  be  minimized  if  the  mesh  is  restricted 
so  that  the  connectivity  is  determined  algebraically.  The  restriction  simply  requires  the 
connectivity  to  be  equivaltiM  tc  a  /cctangular  grid  with  regular  node  numbering.  It  is 
important  to  note  uiat  the  nodal  coordinates  are  not  lestricted  as  long  as  gaps  or  overlaps 
are  not  created.  Voids  can  be  modeled  by  giving  a  rectangular  section  of  elements  zero 
local  stiffness  matrices  and  skipping  the  displacement  solution  at  nodes  that  lie  completely 
within  that  area.  The  connectivity  is  then  calculated  using  the  row  and  column  of  the 
element  However,  if  the  node  numbers  are  determined  solely  by  the  row  and  column 
numbers  the  node  numbers  themselves  do  not  produce  any  new  information.  The  impetus 
for  global  node  uumherin^  is  to  represent  the  assembled  matrices  as  two  chmensiona] 
anays  that  can  lie  inverted  using  stock  solvers.  Since  the  matrices  are  not  inverted,  nodes 
ca'i  be  labeled  using  the-  row  and  column  numbers  by  increasing  the  dimensionality  of  the 
arrays 

U  » 

Where  ^  is  the  displacement  direction  ( 1  =  2  *  y,r) ,  i  is  the  row  number  and  j  is  the 

column  number  This  eliminates  the  operations  necessary  to  convert  to  global  node 
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numbers.  The  connectivity  is  then  easily  calculated  as  the  following. 


Uik,iJ-l)  ,  U(kJ,f) .  Vik.UhD  , 
UikJ*lJ-l)  ,  .  UikJ*lJ*l)  f 


Where  each  term  is  repeated  twice,  first  for  ib «  1  and  second  for  A;  *  2  (see  equation 
(3.10)).  The  elements  are  also  numbered  using  row  and  column  numbers  when  the  unique 
stiffness  and  dancing  rows  are  determined  before  the  first  time  step.  After  that  the 
element  numbers  are  no  longer  needed  using  the  node  by  node  solution. 

Row-Colunui  Multiplication 


To  facilitate  the  row-column  multiplications  in  equation  (4.3),  the  unique  rows  are 
sq)arated  into  two  (3x3)  axrays.  One  array  are  for  the  values  that  are  multiplied  by 
displacements  in  the  x,z  direction  and  the  second  are  for  the  values  multiplied  by  the 
displacements  in  the  y,r  direction.  Equation  (4.3)  becomes 


Fi  •  *  mL,2«4'  *  mL.it',"*' 


Ihc  primed  stiffness  rows  in  equation  (4.9)  are  stored  as  ®  UNKik,l,m,n)  and 

^UNDik,ltm,n)  where  A  indicates  the  displacement  direction,  /and  mare  die  (3x3) 

Hiray  indices  and  n  is  the  unique  row  number.  The  (3x3)  arrays  are  arranged  to 
correspond  to  the  node  connectivity  which  is  alsu  defined  by  a  (3x3)  array  (equation 
(4.8)).  The  row-column  multiplication  can  then  be  accomplished  by  multiplying  values 
in  the  same  array  positions  and  adding  the  result 


Boundary  Conditions 
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Dinchlel  type  boundary  conditions  can  be  applied  simply  by  setting  the  nodal 
displacement  at  time  t  +  At  equal  to  the  specified  displacement.  This  can  be  done  after 
the  solution  at  the  node  is  calculated  or  if  it  is  more  efficient  the  solution  at  the  node  can 
be  skipped.  If  the  Dirichlet  boundary  conditions  are  not  functions  of  time,  the 
displacement  vectors  can  be  initialized  with  the  prescribed  values  and  the  solver  can 
ignoxC  those  nodes.  The  boundary  conditions  are  then  implicitly  applied  at  every  time 
step.  For  the  layered  medium  model,  Dirichlet  boundary  conditions  are  applied  only  at 
the  nodes  along  the  axis  of  symmetry  which  are  restricted  from  moving  in  the  y,r 
direction. 


o 
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Figure  4.3:  Bo>indary  conditions  for  half  symmetry  model. 

Neill  >iann  type  boundary  conditions  can  be  applied  at  the  specified  nodes  by  including  the 
F,  te  m  in  equation  (4.3)  The  force  term  appears  only  at  tlie.  nodes  on  the  top  surface 

undf  c  the  transducer.  This  reduces  the  aniourit.  of  memory  required  to  store  the  array 

and  eliminates  unnecessary  computations  The  remaining  regions  are  traction  free 
boundaries,  sec  Figure  (4.3)  It  is  interesting  to  point  out  that  free  body  motion  is  not 
restricted  in  the  x,z  direction.  This  docs  not  cause  a  problem  because  the  stiffness  matrix 
is  not  inverted.  If  the  forcing  function  has  a  non-zero  net  force,  the  mesh  will  eventually 
gain  a  net  velocity  in  the  same  direction.  The  displacements  will  then  continue  to  grow 
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in  tune.  If  the  forcing  function  has  a  zero  net  force  the  displacements  will  renuiin 
centered  about  the  initial  mesh  location. 


Solution  Regions 

As  the  solver  moves  from  node  to  node  it  needs  to  recognize  whether  the  node  is  on  an 
edge,  on  a  comer  or  in  the  interior.  It  also  needs  to  know  if  boundary  conditions  are 
being  applied  at  that  node.  As  the  solver  visits  each  node  it  determines  the  unique  flag, 
performs  the  row-column  multiplications  and  completes  equation  (4.3)  or  (4.5)  by  adding 
the  remaining  scalar  quantities.  To  accomplish  this  as  efficiently  as  possible,  the  nodes 
are  divided  into  regions  where  the  solution  is  performed  sq)arately.  After  a  number  of 
iterations  the  solver  changes  to  a  difrerent  sectional  configuration. 

The  initial  solver  has  five  sections  as  shown  in  Figure  (4.4).  Section  1  contains  only  the 
comer  node  on  the  top  surface  on  the  axis  of  symmetry.  Section  2  contains  the  nodes  on 
the  top  surface  under  the  transducer.  Section  3  contains  the  nodes  on  the  top  surface  not 
under  die  transducer.  Section  4  contains  die  node  on  the  axis  of  symmetry.  Section  5 
contains  the  interior  nodes,  right  boundary  nodes  and  the  bottom  surface  nodes.  The 
Neumaim  type  boundary  conditions  arc  applied  in  sections  1  and  2.  The  Dirichlet 
boundary  conditions  arc  applied  in  section  4.  The  row-column  multiplications  arc  reduced 
in  sections  1  through  4  because  the  connectivity  is  lowered  from  9  to  6  for  the  edge 
sections  and  from  9  to  4  for  the  comer  node. 
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The 
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solver  also 
has  five 
sections 
(Figure  4.5). 
Section  1 
contains  the 


x,z 

Figure  4.4;  Initial  solver  solution  regions. 


nodes  on  the  top  surface  excluding  the  comer  nodes.  Section  2  contains  the  nodes  along 
the  axis  of  symmetry.  Section  three  contains  all  of  the  interior  nodes.  Sections  4  and  5 
contain  the  tight  hand  side  boundarj  nodes  and  the  bottom  surface  nodes  respectively. 
Similarly  row-column  multiplications  are  reduced  in  sections  1,2,4  and  5. 


Figure  4.5;  Secondary  solver  solution  regions. 


When  computing  the  solution  in  section  5  of  the  first  solver  and  in  section  3  of  the 
second  solver,  the  order  in  which  the  solver  proceeds  from  node  to  node  greatly  affects 
the  execution  tune  for  the  axisymmetiic  formulation.  The  sweep  should  be  carried  out 
by  starting  at  the  upper  left-hand  node  and  moving  down  the  column  and  then  advancing 
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oue  column  to  the  right  This  swcq>  pattern  is  more  efficient  because  the  unique  node 
number  does  not  change  in  a  single  layti  as  the  solution  proceeds  down  the  column.  This 
is  faster  because  the  pointer  to  the  stifi'ness  and  damping  arrays  is  not  recalculated  for 
each  node.  If  a  horizontal  sweep  is  used  epu  time  is  increased  by  30  percent 

Expanding  Solution  Domain 

At  each  time  step  information  travels  across  one  element  For  the  nodes  that  have  not 
been  reached  yet  the  solution  of  equation  (4.3)  consists  of  72  zero  multiplications  and 
72  zero  additions  because  the  displacement  vectors  will  contain  all  zero  entries. 
Obviously  no  new  informaHon  is  gained  by  carrying  out  the  solution  at  those  nodes. 
Therefore  the  solution  will  be  extremely  inefficient  at  early  time  steps  when  a  majority 
of  the  nodes  fall  into  this  group. 

Expanding  the  .solution  domain  at  die  same  rate  as  the  flow  of  information  eliminates  the 
unnecessary  calculations.  The  initial  domain  consists  of  the  nodes  under  the  transducer 
(section  1  and  2  in  the  initial  solver)  and  expands  by  one  row  and  one  column  with  each 
time  step  until  it  reaches  the  boundary  Figure  (4.6)  shows  an  example  of  how  the 
solution  domain  advances. 
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i  c  (•  Figure  4.6:  Example  of  expanding  solution  domain  procedure. 

dnniain  is  most  pronounced  when  the  total  time  of  the  analysis  is  short  As  the  number 

of  iterations  beyond  the  time  when  the  solution  covers  the  entire  domain  incTeases,  the 

realized  time  savings  become  a  smaller  percentage  of  the  total  execution  time.  It  is 

exUoiiieiy  helpful  when  debugging  changes  to  the  program  because  of  the  quick  initial 

response. 

•Solvei  Switching 

Th(  switch  from  the  first  solver  to  the  second  solver  is  determined  by  either  the  number 
of  Iterations  necessary  to  complete  the  solution  domain  expansion  or  the  number  of 

over  which  the  forcing  function  F,  is  applied.  The  greatei  of  ilie  two 

'•  ;nir;'s  Uic  iteration  a!  which  ilie  switch  is  made.  If  thr  expansion  is  the  determining 
>1 .  ihe  forcing  function  wil’  be  sei  to  zero  after  it  is  complete  If  the  opposite  is  true, 
•he  -j  .pansion  will  stop  when  the  outer  boundaries  have  been  reached.  The  solver  will 
not  distinguish  between  the  interior  nodes  and  outer  nodes  when  sections  3, 4  and  5  reach 
tile  edge  of  the  domain  and  will  assume  a  connectivity  of  six,  six  and  nine  respectively. 
The  displacement  vectors  are  dimensioned  to  include  one  extra  row  and  one  extra  column 
of  ‘imaginary’  nodes  on  the  right-hand  and  bottom  sides.  The  solution  is  unaffected 
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because  the  corresponding  stiffness  matrix  con^nents  will  be  zero  and  solutions  are 
never  calculated  at  those  nodes. 

IMsplacement  Vector  Rotation 

The  solution  of  equation  (4.3)  requires  three  separate  displacement  vectors  to  be  stored 
simultaneously,  namely  UO*U*~^  ,  U1  ^D*  and  U2-U**^*.  As  the  time  step 
advances  the  I7‘*^  vector  becomes  D*  and  D*  becomes  Therefore  at  the  end  of 
every  time  step  UO  must  be  set  equal  to  U1  and  then  U1  must  be  set  equal  to  U2.  The 
solution  is  then  repeated  solving  for  the  new  U2.  In  order  to  avoid  this  operation,  die 
solver  is  repeated  three  times  in  succession.  In  the  first  section  UO  =  U1  »  V*  and 
U2^V**^‘.  In  the  second  section  =  U1  and  U2^U\  TTie  third 

section  completes  the  cycle  with  UO  =  f/^  U1  =  U^*^  and  U2  This  process 

shifts  the  displacerr'-  . '  jctors  back  iirtplicitiy. 

If  damping  L  neglected  and  equation  (4.S)  is  used,  only  two  displacement  vectors  need 
to  be  stored.  The  U**^'  displacement  vector  can  be  calculated  directly  into  the 
because  of  the  one  to  one  correspondence  in  equation  (4.5).  In  this  case  the  solver  is 
only  split  into  two  sections.  In  die  first  section  U1  -  W  and  UO  =  to  start  with  UO 
becoming  as  the  solution  proceeds  from  node  to  node.  The  second  section  has 
UO  -  U'  and  U1  =  with  UJ  becoming  U^*^'  to  complete  the  cycle.  For  simple 
models  where  the  number  of  unique  nodes  is  small,  the  displacement  vectors  become  the 
limiting  memory  requirement  For  the  500  by  1000  node  model,  each  displacement 
vector  takes  4  Mbytes  of  memory.  By  rolling  U^*^’  into  the  total  memory  required 
is  reduced  by  nearly  33  percent 
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5  FINITE  ELEMENT  CODE  VERIFICATION 

5.1  Convergence 

Several  factors  influence  solution  convergence.  The  stability  of  the  time  int^iation,  as 
discussed  previously  in  section  4.3,  requires  the  Courant  number  to  be  less  than  or  equal 
to  one.  However,  the  solution  will  not  converge  for  q  >  1.0 .  The  numerical  integration 
of  the  stiffness  matrices,  time  integration  methods,  computer  roimd-off  errors  and 
dispersive  effects  are  possible  causes  of  solution  divergence.  The  inability  of  the  model 
to  represent  the  forcing  function  either  in  time  or  space  will  not  affect  solution 
convergence.  The  solution  will,  however,  not  produce  the  correct  displacements.  The 
mass  matrices  are  not  a  source  of  error  because  they  are  exactly  integrated  by  a  four 
point  Gauss  quadrature. 

Stiffness  Matrix  Numerical  Integration 

The  numerical  integration  of  the  stiffness  matrix  for  the  plane  strain  case  is  exact  for  4 
point  quadrature.  Increasing  the  order  of  integration  will  not  increase  the  accuracy  of 
the  solution.  Single  point  integration  is  not  accq>table  because  it  allows  certain 
combinations  of  displacements  to  exist  in  the  absence  of  nodal  forces.  The  axisymmetric 
formulation  of  the  stiffness  matrix  involves  a  term  which  can  not  be  exactly 
integrated  using  Gauss  quadrature.  However,  for  the  axisymmetric  case  increasing  the 
order  of  numerical  integration  to  9  and  16  point  quadratures  had  no  appreciable  effect 
on  solution  convergence  or  wave  velocity.  One  reason  for  this  is  that  the  wave 
propagations  investigated  were  all  axial  in  nature.  The  axial  waves  are  driven  by  the 
terms  in  the  stiffness  matrix  that  do  not  contain  the  term  and  are  integrated  exactly 
by  the  4  point  quadrature.  For  a  complete  study,  radial  wave  propagation  properties 
must  be  investigated. 
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Hme  Intesration 


The  central  difference  approximation  of  equation  (3.16)  does  not  produce  any  instabilities 
as  long  as  the  Courant  constraint  is  satisfied.  The  backwards  difference  approximation 
in  equation  (3.17),  however,  does  introduce  an  instability  proportional  to  die  magnitude 
of  the  components  of  the  damping  matrix.  The  addition  of  viscous  damping  usually 
increases  the  stability  of  a  numerical  solution  by  damping  out  spurious  modes.  With  the 
backward  difference  formula  the  vdodty  at  time  t  is  estimated  as  the  velocity  att  - 
which  can  cause  the  viscous  force  to  be  larger  than  the  sum  of  the  other  forces  acting  on 
the  node  causing  the  solution  to  diverge  [25].  Decreasing  die  Courant  constraint  will 
decrease  the  time  step  and  reduce  die  effect  of  the  error  of  the  backward  difference 
approximation.  If  the  velocity  was  approximated  with  a  central  difference  formula  the 
viscous  damping  would  indeed  increase  the  stability  of  the  solution. 

Computer  Round-off  Error 

The  code  was  written  in  both  single  precision  (8  digits)  and  double  precision  (16  digits). 
The  solution  convergence  vas  not  affected  by  the  increase  in  accuracy  of  the  double 
precision  representation.  Divergent  solutions  showed  exactly  the  same  A-scan  for  single 
precision  as  for  double  precision.  Double  precision  is  only  necessary  to  retain  accuracy 
during  matrix  inversions.  Because  there  are  no  matrix  inversions  in  the  formulation 
single  preciaon  accuracy  is  acceptable.  However,  because  single  precision  is  more 
sensitive  to  underflow  conditions,  i.e.  when  numbers  become  smaller  than  can  be 
represented  by  the  computer,  equations  (4.3)  and  (4.5)  must  be  modified.  Factoring  out 
the  Ar^  term  from  and  [AT]^  will  nunimize  the  number  of  underflow  checks. 
The  main  advantage  of  using  single  precision  is  a  25  percent  increase  in  processing 
speed. 
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Spatial  and  Temporal  Discretization 


The  degree  of  spatial  and  temporal  discretization  is  determined  by  the  frequency  of  the 
forcing  function.  For  sinusoidal  input,  the  (tegree  of  ^tial  discretization  is  defined  by 
the  number  of  nodes  per  wavelength,  'Hre  degiv  of  temporal  discretization 

is  defined  by  the  number  of  time  steps  per  period,  n,  «  — .  The  degrees  of 
discretization  are  not  independent  due  to  Ae  Courant  stability  constraint.  Substituting 
and  n,  into  equation  (3.20)  yields  the  following  relation 

n,  «  qn,  (5.1) 


which  guarantees  that  the  temporal  discretization  will  always  be  greater  than  the  spatial 
discretization.  The  forcing  function  is  a  raised  cosine  function  given  by 

F  -  [  1  -  cos(— ))cos(«r)  0 i  f  i ^  (5.2) 

3  (•> 

which  has  a  Gaussian  frequency  distribution  centered  at  u,  see  Figure  (S.l).  Figures 
(5.2-4)  compare  displacement  results  for  three  different  values  of  n,  for  the  initial  wave, 
the  first  reflection  and  the  second  reflection,  respectively.  Note  how  the  error  for 
R,  ■  4  does  not  increase  as  it  propagates.  Depending  on  the  degree  of  accuracy 
required,  an  in  the  range  of  6  to  10  is  suggested.  The  numerical  model  used  is 
explained  in  the  next  section. 

5.2  Dispersion 

Dispersion  is  the  variation  of  wave  velocity  with  frequency.  If  the  forcing  function  has 
a  broad  band  frequency  content,  dispersion  will  cause  the  wave  to  change  shape  as  it 
propagates.  An  isotropic,  homogeneous,  lossless  material  is  non-dispersive.  Numerical 
solutions  will  often  introduce  an  artificial  dispersion  creating  an  enor  in  the  solution. 
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Figure  i.Z:  Iziitial  waveform  for  n,  s  4,  7  and  10 
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The  numerical  solution  of  equation  (2.1)  employs  both  qiatial  and  temporal 
discretizations.  Each  of  the  approximations  causes  dispersion.  As  shown  by  Krieg  and 
Key  [27]  and  Belytschko  and  Mullen  [28,29],  choosing  the  sq>propiiate  combination  of 
matrix  and  time  integration  can  mininuze  die  artificial  disperaon.  The  combination 
of  lumped  mass  matrix  and  central  difference  time  integration  for  q  *  1.0  is  a  completely 
nondispersive  solution  along  element  edges.  The  lumped  mass  matrix  decreases  die 
phase  velocity  as  frequency  increases  and  die  central  difference  integration  increases  the 
phase  velocity  as  frequency  increases.  As  q  is  reduced  the  dispersive  effect  of  the  central 
difference  integration  is  reduced  which  shifts  the  net  dispersive  curve  towards  the  lumped 
dispersive  curve.  Because  the  dispo^ion  is  dependent  on  the  element  size  through 
the  Courant  constraint,  a  mesh  with  nonuniform  elements  in  the  same  material  will  have 
nonuniform  dispersive  properties.  For  models  containing  more  than  one  material,  the 
spatial  discretization  must  be  the  same  in  each  material  for  uniform  dispersive 
properties. 

Although  a  completely  nondispersive  solution  is  desirable,  it  is  not  possible  because  the 
solution  diverges  if  extremely  high  frequencies  are  not  attenuated.  In  feet,  for  the 
undamped  case,  solution  convergence  is  completely  determined  by  the  combination  of  the 
highest  unattenuated  ffequencies  and  the  lowest  natural  fluency  at  a  node. 

A  special  model  is  created  to  reveal  die  effect  of  dispersion  on  the  solution.  Starting 
with  the  same  model  as  shown  in  Figure  (4.3)  the  boundary  conditions  on  the  right-hand 
side  are  changed  to  model  another  axis  of  symmetry,  see  Figure  (S.S).  The  forcing 
function  is  applied  to  die  entire  top  surface  with  a  uniform  distribution.  The  raised 
cosine  signal  imparts  a  zero  net  force  on  the  model  and  approximates  the  actual  wave 
pulse  generated  by  a  piezoelectric  transducer  PO].  The  model  dien  becomes  one 
dimensional,  which  eliminates  any  geometric  dispersion.  Figures  (5.641)  show  how  the 
artificial  dispersion  affects  the  wave  shape  as  it  bounces  back  and  forth.  When  q  is  close 
to  1.0  the  solution  is  nearly  non-dispersive  as  predicted  P9];  as  the  q  is  lowered  the 
dispersive  effect  is  increased.  With  decreasing  q,  the  solution  converges  to  die  dispersive 
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ngure  5.6:  Effect  of  ertificial 


dispersion  for  q  «  .995 


Normalized  Displacement  Normalized  Displacement 


Figure  5.7:  Effect  of  artificial  dispersion  for  q  «  .800 


Figure  5.8;  Effect  of  artificial  dispersion  for  q  =  .500 
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e^ect  caused  purely  by  the  mass  lumping. 


5.3  QualitatiTe  Results 

The  wavefronts  produced  by  the  finite  element  code  are  compared  with  Figure  (2.1)  and 
(2.2)  qualitatively  by  plotting  the  displacement  fields  for  a  point  source  and  a  one  quarter 
inch  transducer  at  a  given  time  step.  Surface  plots  are  used  to  di^lay  a  single 
displacement  component.  Although  the  wavefronts  can  be  clearly  identified,  the  plots 
do  not  represent  the  actual  wavefront  shapes  which  contain  only  in-plane  displacements. 
A  7.0  by  7.0  mm  mesh  with  Aluminum  material  properties  is  used  and  displacements  are 
stored  at  1.2  ps  for  the  point  source  and  at  0.8  ps  for  the  transducer. 

Line  and  Point  Source 

Figures  (5.9)  and  (5.10)  show  the  point  force  displacement  results  for  the  plane  strain 
and  axisymmetiic  codes,  respectively.  The  wavefronts  agree  with  those  predicted  in 
Figure  (2.1).  The  longitudinal,  shear,  head  and  Rayleigh  waves  can  all  be  identified. 
The  longitudinal  wave  does  not  fully  ^pear  in  either  plot.  The  x,z  di^lacements  are 
zero  in  the  y,r  direction  and  the  y,r  displacements  are  zero  in  the  x,y  direction  as 
expected,  which  confirms  that  the  displacement  direction  is  parallel  to  the  direction  of 
wave  propagation.  There  is  a  small  bias  in  amplitude  to  the  x,z  direction  which  is 
justified  because  that  is  the  direction  of  the  tq^plied  load.  The  slope  of  the  shear  wave 
in  the  y,r  direction  is  not  zero  at  the  axis  of  symmetry.  This  does  not  create  a 
discontinuity  because  in  the  mirror  image  the  y,r  displacements  must  be  reversed  in  sign 
to  maintain  a  symmetric  response.  There  is  no  sign  reversal  for  the  x,z  displacements. 
The  Rayleigh  wave  wiU  also  have  a  similar  sign  reversal  for  its  y,r  displacements.  The 
Rayleigh  wave  has  the  largest  amplitude  and  does  appear  to  decrease  exponentially  away 
from  the  surface.  The  x,z  and  r,y  displacements  for  the  Rayleigh  wave  are  out  of  phase 
by  90  degrees  which  agrees  with  the  expected  ellipsoidal  motion.  The  large  Rayleigh 
wave  obstructs  the  view  of  the  head  wave  that  originates  where  the  longitudinal  wave 
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Figure  5.9b:  U,  wrface  plot  for  «  Une  ■ouree  ,  /-  5  MHz  i»l2  ftS 


grazes  the  top  surface  and  connects  tangentially  to  the  shear  wave.  The  axisymmetric 
wave  amplitudes  decrease  faster  than  the  plane  strain  wave  amplitudes  as  the  waves 
move  away  from  the  center  because  of  the  difference  in  geometric  dispersion. 

Transducer 

Figures  (S.ll)  and  (5.12)  show  the 
displacement  plots  for  the  one  quarter 
inch  transducer.  Again  the  wavefronts 
agree  well  with  those  predicted  in  Figure 
(2.2).  Most  of  the  energy  from  the 
center  of  the  transducer  goes  directly  into 
the  flat  longitudinal  wave  giving  it  the 
greatest  amplitude.  The  other  waves  are 
smaller  because  they  are  generated  only 
from  the  discontinuity  at  the  edge  of  the 
transducer.  The  waves  caused  by  the  discontinuity  are  similar  but  not  identical  to  those 
caused  by  a  point  force.  The  shear  and  Rayleigh  waves  are  now  inverted  in  the  x,z 
displacements  and  have  the  same  sign  in  the  y,r  displacements.  Figure  (5.13)  shows 
more  clearly  how  the  step  discontinuity  differs  from  the  point  source.  Again,  waves 
moving  away  from  the  center  are  attenuated  faster  in  the  axisymmetric  case  because  of 
the  difrerence  in  geometric  dispersion.  Waves  that  are  travelling  towards  the  center  for 
the  axisymmetric  case  will  actually  gain  amplitude  as  the  spherical  wavefront  collapses. 


Transducer 

fTTfTTT 


Figure  5.13:  Point  source  and  edge  effect 
shear  wave  propagation. 


5.4  Analytic  Results  Comparison 

Analytic  solutions  for  a  line  source  and  a  point  source  at  points  within  the  domain  are 
plotted  against  the  numerical  solutions  as  a  means  of  validating  the  finite  element  code. 
The  analytic  solutions  presented  here  were  calculated  with  the  Cagnaird  de-Hoop 
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Figure  5.11b:  U,  furfKe  plot  for  plane  atnin  witb  a  1/4  inch  transducer,  /*  10  MHz  t  •  0.8  /ar 
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formalism  presented  in  [4]  by  S.  Dai.  A  - r^y.*’ 

raised  cosine  signal  at  10  MHz  is  used  as 
an  input  and  the  points  were  taken  at  radii 
of  R*  1.0,  2.0  and  4.0  mm  each  at 
angles  of  0,  30,  60  and  87  degrees, 

see  Figure  (5.14).  The  results  are  shown  yigure  5.14;  Location  of  points  for  analytic 
in  Figures  (5.15-21).  At  4r-0,  y,r  comparisons. 

direction  displacements  are  not  plotted  because  they  are  set  to  zero  by  the  tilled 
boundary  condition  along  the  axis  of  symmetry.  Reflections  from  the  finite  sized 
numerical  model  are  not  shown  because  the  analytic  solution  assumes  an  infinite  medium. 
All  three  R  locations  are  plotted  on  the  same  graph.  The  R  -  2.5  and  4.0  mm  signals 
have  been  offset  on  the  abscissa  by  1.0  and  2.0,  respectively,  to  separate  the  signals 
completely. 


The  line  source  results  for  the  plane  strain  code  compare  very  well  with  the  analytic 
solutions  for  angles  of  less  than  30  degrees.  As  lir  increases,  the  solutions  no  longer 
agree.  This  is  due  to  the  singularity  of  the  analytic  integration  at  the  x  >  0  surface  which 
distorts  the  head  wave.  Table  (5.1)  lists  the  correct  arrival  time  for  the  longitudinal, 
shear  and  head  waves.  The  bead  wave  is  tangent  to  the  shear  wave  at  4^  *  33.4  degrees, 
therefore  head  wave  timings  are  only  listed  for  60  and  80  degrees.  The  longitudinal  and 
shear  waves  arrive  at  the  same  time  for  all  values  of  if  for  a  given  R.  The  numerical 
solution  gives  the  correct  arrival  times.  The  solutions  also  appear  to  agree  better  once 
the  longitudinal  and  shear  waves  have  separated  at  R  »  2.0  and  4.0  mm.  The  difference 
at  R  •  1.0  mm  might  be  explained  by  a  slight  difference  in  timing  which  causes  the 
waves  to  add  differently  making  the  displacements  appear  to  be  dissimilar. 


The  point  source  numerical  solution  results  are  plotted  in  Figures  (5.22-28).  Although 
the  analytic  data  was  not  available,  the  results  compare  weU  with  those  presented  in  [24]. 
The  analytic  solution  for  a  dirac  delta  function  point  load  is  a  dirac  delta  displacement 
wave,  therefore  a  longitudinal  wave  created  by  convolving  a  dirac  delta  ftmction  into  a 
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Normalized  Displacement  Normalized  Displacement 


Figure  5.15;  Line  source  x-direction  displacement  V'  *  0  /  -  10  MHz 
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Figure  5.16;  Line  source  x-direction  displacement  =  30  /  -  10  Mhz 
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Normalized  Displacement  Normalized  Displacement 
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Figure  5.17;  Line  source  y-direction  displacement  ^  =  30  /  -  10  MHz 


Figure  5.18;  Line  source  x-direction  displacement  =  60  /  -  10  MLz 
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Figure  5.22:  Point  source  z-direction  displacement  V'  =  0  /  -  10  lihz 
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Normalized  Displacement  Normalized  Displacement 


Figure  5.23;  Point  source  z-iJirection  displacement  ^  =  30  /  -  10  MHz 


Figure  5.24;  Point  source  r-direction  displacement  ^  =  30  /  -  10  Mhz 
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Figure  5.25:  Point  source  z-direciion  displacement  -  60  f  -  10  MHz 
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Figure  5.26;  Point  source  r-direction  displacement  ^  =  60  /  -  10  Mhz 
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Normalized  Displacement  Normalized  Displacement 


Figure  5.27:  Point  source  z-direction  displacement  V*  =  B7  /  -  10  MHr 


Figure  5.2B.  Point  source  r-direction  displacement  ^  -  B7  /  -  10  Mhz 
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Table  5.1:  Longitudioal,  shear  and  head  wave  arrival  tiines  in  /is. 


1  ndiiis 

Long. 

Head  Wave  | 

Shear 

60* 

87*  1 

1.0 

.1597 

.2594 

.1722  1 

2.5 

.3993 

.7251 

.6484 

.4304  1 

4.0 

.6389 

1.1601 

1.0374 

.6887  1 

raised  cosine  produces  a  raised  cosine  displacement.  The  longitudinal  wave  for  the 
numerical  solution  is  a  raised  cosine  as  expected.  The  wave  timings  also  agree  well  with 
those  listed  in  Table  (5.1). 

5.5  Experimental  Results  Comparison 


A  total  debond  is  modeled 
experimentally  using  a  single 
aluminum  plate  supported  at  the 
perimeter.  The  transducer  is  placed 
in  the  center  of  the  plate  to  avoid 
contaminating  the  A-scan  results  with 
reflections  from  the  plate  edges.  One 
quarter  inch  diameter  transducers  with 
center  frequencies  of  S  and  10  MHz 


Figure  5.29:  Single  aluminum  plate 
experimental  configuration,  (not  to  scale) 


are  used.  Figure  (5.29)  shows  the  physical  dimensions  of  the  experimental  configuration. 
The  material  properties  of  the  aluminum  are  Cj> 6261.0  m/5,  3448.2  m/5  and 

p  *  2842.0  kgim} .  The  numerical  A-scans  are  calculated  by  a  weighted  average  of  the 
displacements  at  the  nodes  under  the  transducer  for  each  time  step.  The  weighting  is  the 
same  as  the  distribution  of  the  uniform  load  to  the  nodes. 


The  experimental  and  numerical  results  are  compared  in  Figure  (5.30).  The  experimental 
data  is  shifted  so  that  the  first  longitudinal  reflection  overlays  with  the  numerical 
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Figure  S.30a:  N\imerical  and  experimental  A-seans  for  debond.  /  -  5  IfHz 


Time  (microseconds) 


Figure  5.30b:  Numerical  and  experOTcntal  A-scans  for  debond,  /  =  10  kHz 
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prediction.  The  axisymmetiic  code  without  material  damping  is  used  for  the  numerical 
solution.  The  experimental  and  numerical  results  are  in  excellent  agreement  even  without 
any  damping  effects.  The  large  waveforms  are  the  first  and  second  longitudinal  wave 
reflections  from  the  bottom  surface.  The  shallow  waveforms  between  the  longitudinal 
reflections  are  mode  converted  shear  waves  created  when  die  longitudinal  wave  reflects 
from  the  bottom  surface.  Inspection  of  the  experimental  results  suggests  that  the  input 
wave  is  dose  to  the  raised  cosine  but  not  exactly  the  same. 


IMMi 


— 

Aluminum 

1  ^  AooutSeM 

r 

1 

Jmm 

Aluminum 

1 

■ 

■ 

-y/ 


A  schematic  of  the  layered  material 
experimental  configuration  using  two 
aluminum  plates  separated  by  a  thin 
layer  of  acoustic  gel  is  shown  in 
Figure  (5.31).  The  material  ^ 

properties  of  the  gel  arc 
C^•  1490  m/s  and  1080.0  kglm^. 

The  shear  wave  velodty  for  the  gd  is  figure  5.31;  Three  layered  medium 
veiy  tow  and  is  estimated  for  the  configmation.  (not  to  scale) 

numerical  solution  as  10  percent  of  C^.  The  bond  thickness  is  controlled  by  a  small 
diameter  wire  placed  around  the  perimeter  of  the  plates.  Physical  dimensions  are  also 
shown  in  Figure  (5.31). 


Figure  (5.32)  compares  the  experimental  results  to  numerical  solutions.  The  first  wave 
is  the  longitudinal  wave  reflection  from  the  bottom  of  the  top  plate.  The  following 
smaller  wave  is  the  longitudinal  wave  reflected  from  the  top  surface  of  the  lower 
aluminum  plate.  The  next  two  waves  are  longitudinal  waves  that  have  reflected  back  and 
forth  within  the  bond  layer  multiple  times.  The  phase  difference  in  the  bond  layer 
reflections  could  be  caused  by  the  lack  of  damping  in  the  numerical  model. 
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Figure  5.32a:  Numerical  and  experimental  A-scans  for  .466  nun  bond  /  «  5  l£Hz 


Figure  6.32b;  Numerical  and  experimental  A-scans  for  .466  mm  bond  f  -  10  MHz 
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5.6  Attenuation  Example 


The  damping  effect  on  wave  propagation  is  illustrated  by  plotting  die  di^lacements  at 
successive  time  intervals  with  the  same  scale.  The  z  direction  displacements  are  plotted 
for  a  undamped  and  a  damped  case  at  f  *  1.0»  2.0,  3.0,  and  4.0  fts  in  Figures  (5.33-36). 
The  damping  is  assumed  to  be  isotrqiic  with  values  of  •  111.45  NS/m^  and 
Dp,, » 25.07  NS/m^  [11]  so  that  the  damj^g  matrix  is  proportional  to  the  stif&iess 
matrix.  This  situation  is  equivalent  to  Rayleigh  damping  with  a  «  0  [16]. 

The  attenuation  of  the  longitudinal  wavefront  is  evident  after  only  one  reflection.  At 
r  >  4.0  iis  the  longitudinal  wave  has  been  reflected  by  the  top  surface  for  the  first  time. 
Introducing  viscous  damping  through  equation  (2.2)  results  in  damping  that  is 
proportional  to  the  square  of  the  frequency,  which  results  in  a  dispersive  solution. 
Therefore,  in  addition  to  decreasing  the  amplitude  of  the  waves,  the  damping  term 
produces  a  change  in  the  waveform.  The  actual  damping  in  a  polycrystaUine  solid  such 
as  aluminum  has  a  more  complex  attenuation  that  becomes  proportional  to  the  frequency 
to  the  fourth  power  for  very  high  frequencies  [31]. 

5.7  Anisotropy  Example 

The  anisotropic  capability  of  the  code  is  tested  using  the  transversely  isotropic  [Q 
matrix  for  Cobalt.  The  components  of  [C]  are 

C„j,  -  35.81x10"  ^yw* 

Cjjjj  »  10.27x10"  N/m^ 

-  1630x10“  Nim^  (5  3) 

C2222  -  30.70x10"  N/m^ 

^1112  “  7310x10"  ^ym’ 
p  »  8900.0  kg/m^ 

Figure  (5.37)  shows  the  z  and  r  direction  displacements  respectively  at  /  -  1.8  ps  for 
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Figure  5.37a:  U,  surfue  plot  for  transversely  isotropic  Cobalt,  /*  S  MHz  t  ~  1.8  iXS 


Figure  5.37b:  U,  surface  plot  for  transversely  isotropic  Cobalt,  /  =  5  MHz  t  =  1.8  /XT 


66 


a  point  source  with  a  S  MHz  signal.  C^.,  in  the  Courant  constraint  is  .  The 

results  are  the  same  as  those  published  [11].  The  anisotropic  analy^  introduces  more 
error  into  the  solution  because  the  wave  velocity  is  now  direction  dependent.  As  a 
result,  the  Courant  constraint  varies  with  direction  introducing  even  more  complicated 
dispersive  mesh  properties. 
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INTRODUCTION 

Adhesive  bonding  has  found  extensive  application  in  the  aircraft  and  defense 
industries  where  the  failure  of  a  bond  in  any  of  the  critical  load-bearing 
components,  for  example  the  rotor-blade  of  an  helicopter,  can  bring  about  a 
catastrophic  failure.  Nondestructive  evaluation  of  adhesively  bonded  structures 
attempt  to  assess  the  key  factors  of  bond  strength  and  quality.  The  bond  strength 
[1]  is  primarily  determined  by  the  thickness  of  the  bondline,  as  this  greatly  affects 
the  stored  energy  in  the  bond.  Three  factors,  if  determined,  provide  a  good  measure 
of  bond  quality.  They  are  bond  thickness,  contact  angle  of  adhesive  to  substrate, 
and  substrate  surface-free  energy. 

The  objective  this  paper  sets  out  to  achieve  is  a  merjis  of  accurately 
determining  the  thickness  of  a  bond  layer.  Digital  signal  processing  is  used  to 
analyze  the  reSected/transmitted  data.  Both  the  time  domain  and  frequency 
domain  approaches  are  investigated.  In  the  frequency  domain  the  Chirp-Z 
transform  [2], [5]  is  employed  to  perform  the  spectral  analysis  of  the  received  signal. 
Based  on  the  resonance  effect  of  the  bondline  the  dips  in  the  observed  transducer 
frequency  spectra  can  be  correlated  to  the  thickness  of  the  bond.  This  topic  is 
further  discussed  in  references  [6]  and  [7]. 

As  an  alternative,  time  domsiin  analysis  of  the  reflected  signal  consists  of 
modeling  the  bonded  structure  as  an  adaptive  Alter  [3],  with  the  taps  being 
representative  of  the  adhesive  layer.  This  technique  is  used  to  provide  a  deconvolved 
bondline  response.  Adaptive  filters  have  found  wide  use  in  adaptive  modeling  and 
system  identiflcation  in  which  the  unknown  system  is  described  by  its  input-output 
behavior.  The  adaptive  filter  then  tries  to  emulate  the  system  response  due  to  a 
known  input.  In  its  application  to  deconvolution,  the  effect  of  undesired  medium 
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Figure  1:  Generic  configuration  used  for  the  experiments 

a.  TVansducer  in  pulse  echo  mode 

b.  Substrate  layers 

c.  Copper  wire  spacers 

d.  Adhesive  layer  (Gel) 

influence  on  the  signal  is  sought  to  be  negated.  Here  the  medium  which  affects  the 
bond-line  response  can  be  considered  as  the  cumulative  effect  of  the  substrates,  the 
transducer  itself  and  the  associated  hardware  that  make  up  the  entire  experimental 
arrangement. 

EXPERIMENTAL  ARRANGEMENT 

The  experimental  arrangement  for  this  research  required  the  assembly  of 
bonded  structures  of  predetermined  thicknesses,  as  shown  in  Figure  1.  The 
substrate  layers  &,  are  optical  glass  flats  of  precisely  measured  thickness.  They  are 
approximately  isotropic  in  nature  with  known  material  constants.  Knowledge  of 
material  parameters,  like  the  density  and  the  longitudinal  velocity  of  sound  in  each 
of  the  three  layers  is  essential  in  determining  the  thickness  of  the  adhesive  layer. 
The  adhesive  layer  d,  is  simulated  by  a  viscous  gel  of  known  material  properties.  To 
corroborate  the  predictions  obtained,  copper  wires  c,  of  known  diameters  were 
placed  as  spacers  in  between  the  substrates.  The  contact  transducer  (0.25  inch 
aperture)  placement  a,  is  generally  in  the  Pulse  Echo  (P.E)  mode  ,  with  a  chosen 
range  of  center  frequencies  between  5-20  MHz. 


Figure  2:  Experimental  Arrangement 
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ADAPTIVE  FILTERING  APPROACH 

Adaptive  filters  are  based  mainly  on  three  diflerent  approaches:  the  Wiener 
Filter  Theory ^  the  Kalman  Filter  Theory ,  and  the  classical  method  of  least  squares. 
Since  the  Wiener  and  Kalman  filters  involve  a  stochastic  formulation  the 
deterministic  least  squares  approach  is  chosen.  This  method  requires  minimizing  an 
index  of  performance  related  to  the  error.  The  error  is  defined  as  the  difference 
between  the  desired  response  and  the  actual  filter  output.  Three  different  classes  of 
algorithms  form  the  basis  of  the  method  of  least  squares.  They  are  recursive  least 
squares,  least-squares  lattice  algorithm,  and  the  QR  decomposition  least  squares 
algorithm.  In  this  paper  the  recursw  least-squares  algorithm  forms  the  basis  of  the 
deconvolution  technique. 

RECURSIVE  LEAST  SQUARES  ALGORITHM 

The  structural  basis  of  the  recursive  least  squares  (RLS)  algorithm  is  a 
transversal  filter  illustrated  in  Figure  3,  where  v(n)  represents  the  input  to  the 
filter,  d(n)  represents  the  desirsd  response,  d{n)  is  the  output  of  the  filter,  and  e(n) 
is  referred  to  as  the  estimation  error.  The  RLS  algorithm  starts  by  defining  a 
correlation  matrix  $(n)  such  that 

#(n)  =  ^ A""*u(i)u^(i)  where  0<A<1  (1) 

«=i 

here  A""*  is  an  exponential  weighting  factor,  or  memory  factor.  The  deterministic 
normal  equation  is  defined  as 

♦(n)w(n)  =  ©(n)  (2) 

with  w(nl  being  the  optimum  value  of  the  tap- weight  vector  for  which  the  index  of 
performance  attains  its  minimum  value.  In  equation  (2)  ©(n)  is  the  deterministic 
cro^'i-correlation  vector  defined  as 

®(")  =  (3) 

i=l 
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The  index  of  performance  £,  sought  to  be  minimized  in  this  case  is 


«(")  = 

(4) 

»=i 

with 

e(t)  =  d(0  -  w^(n)u(i) 

(5) 

Rephrasing  equation  (2)  we  can  get  an  expression  for  the  optimal  tap-weight 
vector  as 


wr(Ti)  =  $“'(n)©(n) 


(6) 


The  above  equation  involves  the  inversion  of  an  M-hy-M  matrix  which  is 
accomplished  by  employing  the  Woodbury  [3]  matrix  inversion  lemma. 

APPLICATION  TO  BONDED  STRUCTURES 

Figure  4a  illustrates  the  P.E  test  of  a  triple  layered  medium.  The  reflected  signal 
d{t)  is  interpreted  to  be  the  convolution  of  the  signal  n(t)  tahen  from  a  single  layer 
of  substrate  material  as  shown  in  Figure  4b,  and  the  adhesive  layer,  represented  as  a 
transfer  function  h{t)  of  Figure  4c.  The  signal  u(t)  can  be  considered  to  have  all  the 
external  influences  sought  to  be  negated  ixom  the  bondline  response.  Referring  to 
Figure  3,  it  can  be  seen  that  here  the  input  to  the  adaptive  Alter  would  be  u(t)  and 
the  desired  response  to  which  u(t)  is  adapted,  is  d{t).  On  convergence  the  Alter  taps 
are  then  representative  of  the  bondline  transfer  function  h(t). 


Figure  4:  a.Signal  from  bonded  structure  d(t)  =  u(t)  *  h(t) 

b.  Signal  from  single  layer  of  substrate  material  u(t) 

c.  Bondline  transfer  function  h(t) 
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RESULTS 


The  case  study  presented  here  deconvolves  the  bondline  transfer  function  of  a 
simulated  bond  of  a  measured  thickness  of  117.6  fim.  The  adaptive  filter  set  up  has 
256  taps.  The  input  signal  u(t)  consists  of  512  samples,  and  the  desired  response 
d(t)  also  has  512  samples,  v(t)  is  delayed  by  256  samples  to  improve  performance  of 
the  filter  [4]. 


Figure  5;  Process  of  Adaptation 
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Tap  Weight  (hJnD 


Figure  5  describes  the  details  of  the  process  of  adaptation.  In  this  figure  e(n) 
shows  the  error  signals  at  each  iteration  which  describe  the  convergence  of  the  filter. 
When  the  a  priori  error  and  the  a  posteriori  error  have  the  same  value,  the  filter 
has  converged  to  its  optimum  tap  weight  vector.  The  signal  d(n)  describes  the 
predicted  output  which  closely  matches  d(n).  The  filter  taps  on  convergence 
represent  the  adhesive  layer  transfer  function  depicted  in  Figure  6a.  This  transfer 
function  is  then  convolved  with  an  idealized  transducer  signal  having  a  Gaussian 
frequency  response  centered  at  10  MHz,  shown  in  Figure  6b. 

The  result  of  the  convolution  is  shown  in  Figure  7a  in  the  time  domain.  Figure 
7b  illustrates  the  frequency  domain  response  of  this  signal.  When  the  frequency 
spectrum  of  the  experimentally  obtained  signal  d(n)  is  overlud,  it  can  be  seen  that 
location  of  the  nulls  in  either  case  are  in  close  alignment. 


so  100  ISO  200  2S0  O  .05  .1  .15  .2  .25 

Tap  Number  Time  (micro  secs) 


Figure  6:  a.  Filter  taps  representing  h{n) 

b.  Idealized  10  MHz  transducer  signal  RCo3{n) 
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Figure  7:  &.Time  domain  response  of  RCo3{n)  *  h{n) 

b.Frequency  domain  response  of  RCos(n)  *  h{n) 
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Finally  in  order  to  validate  the  deconvolution  model  presented  here,  v(n)  and 
h(n)  are  convolved  and  the  results  compared  with  d{n).  Figure  8a  is  the  time 
domain  response  of  this  convolution  which  compares  favorably  with  d(n)  shown  in 
Figure  5.  The  comparison  in  the  frequency  domain  of  the  signals  tt(n)  »  h{n)  and 
d(n)  is  shown  in  Figure  8b.  As  can  be  seen  there  is  an  almost  one  to  one 
correspondence. 


Figure  8  a.Time  domain  response  of  u(n)  *  h(n) 


b.Frequency  domain  response  of  u(n)  *  h{n) 
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CONCLUSION 

This  paper  discusses  the  precise  thickness  meuurement  of  thin  adhesive  layers 
by  employing  a  method  of  deconvolving  bondline  response.  For  the  thickness 
measurement  a  physical  model  was  developed,  to  which  a  frequency  domain  analysis 
can  be  applied.  Successful  measurement  of  thin  layers  based  on  this  frequency 
domain  model  has  been  demonstrated  [6], [7].  To  deconvolve  the  adhesive  layer 
response  a  transfer  function  model  is  developed  based  on  adaptive  filters.  The 
deconvolved  bondline  transfer  function  is  obtained,  and  compared  to  experimentally 
recorded  signsis.  These  comparisons  show  exceUent  agreement. 

The  application  of  this  technique  to  realistic  bonding  configurations  is  one  of 
the  future  goals.  Identification  of  broad  band  signals  to  replace  the  simulated 
transducer  signal,  would  provide  a  better  frequency  evaluation  of  the  bondline 
transfer  function.  Future  work  may  also  entail  applying  the  transfer  function  model 
to  pattern  recognition  schemes. 
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APPENDIX  A 


The  elastic  constitutive  matrices  for  isotropic  materials  in  terms  of  Lam6’s  constants 
are  shown  in  the  index  contracted  fcrms  of  the  stress  strain  relations  for  plane  strain 
in  equation  (A.l)  and  for  axisymmetry  in  equatioi  (A.2). 
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